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A.V. BALAKRISHNANf 


Abstract 

The temporal spectral density of the log-amplitude scintillation of a laser 
beam wave due to a spatially dependent vector-valued crosswind (deter- 
ministic as well as random) is evaluated. The path weighting functions for 
normalized spectral moments are derived, and offer a potential new tech- 
nique for estimating the wind velocity profile. The Tatarskii-Klyatskin 
stochastic propagation equation for the Markov turbulence model is used 
with the solution approximated by the Rytov method. The Taylor “frozen- 
in” hypothesis is assumed for the dependence of the refractive index on 
the wind velocity, and the Kolmogorov spectral density is used for the 
refractive index field. 


f Research supported in part under grant NCC 2-374, NASA. 
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1. Introduction 


This paper is in two parts — this first part deals exclusively with theory. The 
second part will be devoted to computational/experimental results. 

The study of the log-amplitude scintillation spectrum of a coherent beam in wind 
turbulence was initiated largely in the late-1960 to early-1970 period spurred by 
the application to “remote sensing” — see [7-15]. We are particularly interested in 
the case where the wind velocity varies along the propagation path — as in “wind- 
shear,” of importance in flight systems. Most of the work involves calculating the 
space-time correlation function for a spherical wave or a plane wave. The “frozen- 
in” Taylor hypothesis is used to account for the (cross) wind, and the Kolmogorov 
spectral density is used for the refractive index field. Lee and Harp [8] calculate the 
correlation function for the spherical wave while Ishimaru [2] calculates the spectral 
density for a constant wind velocity, in the plane wave case. In this paper we calculate 
the temporal spectral density for an arbitrary spatially varying vector wind velocity 
for the general beam wave case. 

To estimate the wind velocity profile Barakat and Buder [14], following the early 
work of Lawrence et al. [9], use the “slope method." The space-time correlation 
function for the spherical wave is shown to be stationary in time and space and 
further the time derivative at time zero is a linear functional of the magnitude of 
the velocity. The corresponding “path weighting function” is calculated, resulting 
in a linear integral equation for the wind velocity (magnitude). In this paper we 
show that the normalized moments of the scintillation spectral density yield a similar 
relationship and calculate the corresponding path weighting functions. The spectral 
moment of order two is of particular importance because of its relation to the zero- 
crossing rate (measurable real time) given by Rice’s theorem [18]. We have thus a 
potential real time alternate technique for estimating wind velocity. 

We also consider a random wind model in which the wind velocity is assumed to be 
a Gaussian 2D random field with given mean and covariance matrix. We calculate the 
spectral density as well as the spectral moments of order two and the corresponding 
path weighting function. 

We begin with a brief review in Section 2. The necessary propagation theory is 
well described in references [1,2,3] and is reviewed further in [4]. We use the Markov 
turbulence model leading to the random parabolic Tatarkskii-Klyatskin propagation 
partial differential equation. We use the generally accepted approximation to the 
solution due to Rytov, and calculate the slow-varying component of the logarithm 
of the amplitude assuming the Taylor hypothesis to account for the wind. Section 3 
contains the spectral density and spectral moment calculations for the deterministic 
case and Section 4 for the random wind model. 
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2. Review of Propagation Theory 

The starting point of all propagation theory (as in [1,2,3]) is the Helmholtz equa- 
tion for the scalar electric field over R 3 for each component E of the (complex) electric 
field (in the usual notation): 

V 2 £ + k 2 n 2 E = 0, n = c 2 pe 


where k is the wave number: 


k = 




wavelength 

and n( ) is the refractive index. Next we choose a “propagation” direction: say, 
X-axis, and write 

r ~ix + jy T kz 
P = jy + kz 

where i, j, k are orthonormal unit vectors in 3-space. Then we let 

E(x,p) = v(x, p)e lkx 


and substituting in the Helmholtz equation we obtain: 

d 2 v dv 9 . 9 

— - T 2 ik— + k 2 (n 2 — l)u + Av -- 0 
ox * .ox 


where 


Now we can write: 


Av = 


' d 2 d 2 \ 
K dy 2 + dz 2 ) V ' 


d 2 v , dv 

+ 2 ' k di = 


We assume that v is “slow- varying” in x, so that 



+ 2ikv 


) 


dv 

1 - 2 ikv ~ 2 ikv, 

dx 

or that the term can be neglected. Then we have the “parabolic” approximation 
[1,3,4,17]: 

J 

2 ik — T Av + k 2 (n 2 — l)v = 0 (2.1) 

ox 

which we may treat as an “initial value” problem in x > 0. 
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If n = 1, (2.1) simplifies to: 


dv 

2 ik — — h Av = 0 
ox 


( 2 . 2 ) 


or 


— - — Au 
dx 2 ik 



x > 0. 


(2.3) 


We can treat this as in [4] as an abstract Cauchy problem for v(x, ■) over the Hilbert 
space Ti = 2/ 3 (Jt a ): 


where 



is the infinitesimal generator with 


V(A) = n 2 . 


The solution to the initial value problem is then given by 

u(x,-) = S(x)v(0, •) 

where S(-) is the semigroup generated by A. As in [4], we can use Fourier transforms 
with _?"(•) denoting Fourier transform for each x: 

?(v(x)i A) = j e MA.Pi v (x ) p )|dp|, A € R 2 

By slight abuse of notation, we shall write p 2 for |p| 2 ; A 2 for |A| 2 . Since 

T(Av,\) = -47r 2 A 2 ^(u; A) 


we have 

JF(v(x);A) = ^(S(x)u(O, •)) = A). 

Hence 

u(x,p) = f G{x,p - p')v(0, p)\dp\ 

J F? 

where the Fourier transform of G(x, p), for each x > 0 is given by: 

e -4n 2 \ 2 (ix/2k) ( 2 . 4 ) 
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Hence 


and 


-1 k 


1 ikp 2 


G(l ' p) = i7||) exp T Tx p ' = 2^ exp 2 x 


v{x,p) = 


lik\p-p'\ 2 , n ^ , t 

. , exp- v(0,p)|dp|. 

27 TlX Jf& 2 X 


L 


(2.5) 


For the special case where u(0, •) is a 6-function at the origin in R 2 , we obtain 

k e (ikx+(kip?)/2x) 


E{x,p) = 


2i \ix 


Since constant multiphers do not matter, we may take 


E(x,„) = —e 


1 1 *<hV/2*)) 


( 2 . 6 ) 


which is the result obtained by Ishimaru [2, p. 377] although by somewhat different 
arguments. 


Beam Wave 

For the beam wave case, the initial condition is (as in [2]) given by: 

—p 2 ak 


E(0,p) = exp 


where 


a = 


7 rW 2 


+ 


ib 


W — beam size ; Rq = focal length 


Ishimaru [2] goes through a similar procedure as in deriving (2.6) to obtain E(x, •) 
but here we again get this directly from our parabolic approximation (2.1). Thus 

E[x,p] = e ikx v(x,p) 

and writing v(x) = v(x, •), 

v(x) = S(x MO); v(0, •) = E(0, •)• 

Using Fourier transforms to evaluate (2.5), we have 


o 



Hence 


W * );A) = (^l) eXP ' 47r2 (^ 

= ( S ) exp - 4 * 2A2 ( 


■= H -* v ( s ) 


<r 2 + rt 


= (|)exp-4^A 2 (i±i 


v{x,p) 


— if otkp 1 


exp — - 

1 + ixa 2 V 1 4* ixa 


OUU y , v / O V 

, . f e lkx \ -if akp 2 \ 

E{x ’ p) = (T7i^] exp T(rn^j 

as in [2, p. 381]. 

From (2.7) it follows as in [2, p. 381] that the beam size at x is given by 

[7 x \ 2 4x 2 

w fv~Rj + i fc ¥' 


Turbulence Equation 

Consider now the case n ^ 1, we have, rewriting (2.1) 

dv Av — k 2 (n 2 — l)u 

dx 2 ik 2ik 

dv i . ik 2 „ 

dx 2k 2 V ’ 


If we take 


n = 1 + n, 


n 2 - 1 = (1 + n) 2 — 1 = n 2 4- 2n 
omitting the n 2 as small compared to 2 n, we have: 


n 2 — 1 = 2n 
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— = —Av + iknv (2.8) 

dx 2k 

where nv denotes the function n(x, p)v(x , p ). 

Exact solution of (2.8) not being possible, we invoke the two well-known approx- 
imation techniques [1,2,3]: the Born approximation and the Rytov approximation. 

Born Approximation 

Rewriting (2.8) as an integral equation: 

v(x) = S(zMO) + ik f 0 S{x - a)(h(a)v{a)) da 

and expressing the solution in a Volterra series expansion (see [4]), and truncating 
the series at the linear term yields the Born approximation [1,2,3]: 

u(x) = S(x)u(0) + ik f S(x - (7)[n(ff)5(o-)v(0)] da 

or, 

v(x, p) = Eq{x , p) + ik J Q J r2 g ( x - a >p~ p')' p') E o( a > p') \ d p '\ > P e r2 - 

Using the next approximation 

log(l + 2 ) = 2 for |x | small, 
we can reduce this further to yield 

x(x,p) = log|u(x,p)| = log \Eo(x,p)\ Re'll) b(x,p) 

where 

ip B {x,p) = „ ) k " T [ [ G(x-a,p- p')h(a,p)E 0 {a } p) \dp\. (2.9) 

Eo(x l p)Jo Jr 2 

We now specialize to two cases: the plane wave and the beam wave. 
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Case 1. Plane Wave 


Eo(x,p) = e 


ikx 


and (2.9) yields 
4>b(x,p) = 


2 r x p ik(e- x ) 


-/ 

27r Jo 


, r i ik\p- p i 2 ^ u /| 

do exp — — — n(a,p)dp| 

(x - a) Jr 2 2 (a - x) 


(2.13a) 


Case 2. Beam Wave 

Here , _ 

e lkx —1 akp 2 

£^o(a:, P = 7 — — exp — — — — 
1 4- iax 2 1 4- iax 

Hence (2.9) yields 

, x -hct ( 1 akp 2 \ f x do 

<, B (x,p) - ifc(l + iax)e J • jf Q — 


A: 

27T i 


f f 1 P - P \ /\ 

/ ( ex P r P T 

Jr 2 y 2 x - a J 1 




— 1 a/c/r 


— — exp — — — — | dp \ 
+ taa 2 1 + tax 


1 ik\p — p /|2 
x — a 

which, as in (2.14), can be expressed 

k 2 f x e 8fc(<T-l) t ik \p' — 7 (p)p | 2 ,, , ,, ^ 

— / 7 77 — : / , exp ~n — 7T7 7 n(<r, p ) |dp |. (2.14a) 

2tt Jo (x — a)'y(a) Jr 2 2 y(a)(x — a) 

Next we consider the Rytov approximation [1,2,3]. 


Rytov Approximation 

We begin by expressing the electric field as: 

E — Eoe^; ^ (0) = 0; E(0) = E o (0) 

in the Helmholtz equation, where Eq is the solution of (2.1) for n = 0, or equivalently 
of (2.3). Then 

(V 2 +n 2 k 2 ){E 0 e' i ’) = 0 
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yields 
We now let 
Then 


(V 2 + k 2 ){E 0 i>) = k 2 {n 2 -\)E 0 . 


( 2 . 10 ) 


E 0 ip = ve ikx . 


' d 2 v 


(V 2 + * V“* = ( + Av + 2.Jr - * 2 ^ 


Akx 


where we invoke “the parabolic approximation” and neglect the term. Then (2.10) 


yields 


2ik^~~ + Av = —k 2 (n 2 - 1 )(E 0 e ikx ) 

OX 


or 


where 


since 


dv iAv 


+ iknEoe 


—ikx 


dx 2k 

v(0) = Eo( O)V'(O) = 0 


(2.11) 


tp{ 0) = 0. 

Hence, solving (2.11) we have the so-called “Rytov First Iteration Solution,” [2], 

p ikx ik r* 

= 4-, —J S(x - a)(h(a, .)) e - to do (2.12) 

B 0 (x, p ) Jo 


where 

n(a, -)E 0 {a, •) 

represents the function 

n(a, p)Eo(o, p), p e R 2 ■ 

We specialize next to two cases. 


Case 1. Plane Wave 
Here we take 

E 0 (x) = e ikx 


and hence 

n(a,p)E 0 {a l p) = n(o,p)e lka 
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and ( 2 . 12 ) yields 


or, 


V»(x.p) = ( ik ) j da J R 


Tp(x) = ik f S(x - o)(n(o, •)) do 
Jo 

1 ik\p - p '\ 2 


ex Po 


i 2 

27T 


2 x — a 
/|2 


n(a,p') |dp'| 


/r2 27u(x - u) 

[ — - — f — - — exp ^ — — n(cr, p') \dp'\. (2.13) 

Jo x — a Jr 2 x — a L x — a 

This agrees with [2, (17-27a)] but we obtain it by the parabolic approximation. 


Case 2. Beam Wave 
For the beam wave case 


E( 0 ,p) = £o( 0 ,p) = exp 


—p 2 ak 


and as we have seen (cf. 2.7)) 


Eo{x,p) = 


Akx 


■ exp 


— 1 akp 2 


1 + ixa 2 1 + ixa 


Hence 


/ 1 4- ixa 

t(x,p) = {ik) ( — 3 — exp 


1 akp 2 \ 

2 1 + ixa J 


k 

27 ri 


[ x do f ( 1 ifcp-p 2 \ 1 -1 akp 2 , 

/ / I ex P r — -)n{o,p )—-. — exp— : — | dp \ 

Jo x — a Jr? \ 2 x - o J \ to a 21 


+ toa 


which as Ishimaru [2] has shown can be simplified to 

k 2 f ik |p' - 7 (tf)p | 2 


rp(x,p) = f 

Jo 


do 


2 n')(o)(x — cr) 


f exp l — — 7 - n(o,p') \dp'\ (2.14) 

Jr? 2 7 (<r)(x - a) 
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where 


y( a ) = YT~ = 7*(<0 + * 7 /(*) 
1 + tax 


7*(<0 = 


7/( <J ) 


: ~ g j— < 0 

(i't) + a R x2 


2 

QR " kW 2 ' 

As is known (see [2, p. 349]) the Rytov approximation is generally accepted as 
superior to the Born approximation, and so we shall use only the Rytov solution in 
the rest of the paper. 


Markov Turbulence Model 

To proceed further we need to specify the stationary covariance function of the 
refractive index field n(-) in the propagation equation in (2.8). Here we follow the 
Tatarskii-Klyatskin theory [1] and invoke the “Markov turbulence model” — that it is 
“delta-correlated” along the propagation direction. Denoting the covariance function 
by Ra(a , p) we have: 

R*(a,p) = CiHe)R(p) (2.15) 

where C 2 in the usual notation denotes the turbulence intensity. To determine the 
function R(-) we again follow [1,2] and require that 

/ OO _ yoo 

R n (a,p)da = / R z (a,p)do 

-OO J -oo 

where R z (a, p) is the covariance corresponding to an isotropic random field, with the 
Kolmogorov spectral density. This yields 

R(p) = jf i e 2l ‘ |A ' ,| Q(|A|) |dA| 

where 

Q(A) = -TT 76 exp- 47 T 2 ^A 2 . ( 2 . 16 ) 

(* + ^ 2 ) 
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It is customary to allow the turbulence intensity to depend on the propagation 
direction so that n(-) is now a “locally stationary process.” However usually little 
is known about this dependence. We shall avoid this uncertainty by restricting con- 
sideration, if necessary, to propagation path intervals small enough so that it may 
be assumed to be constant. Moreover the quantities of interest to us will be “nor- 
malized” as in [9] so that we do not need to know the precise value either. In what 
follows we shall omit this as a multiplicative constant. 


3. Response to Cross Wind 

We are now ready to study the response to wind. Let w(<r), 0 < o < x, denote 
the “crosswind,” i.e., projection of the wind velocity in the plane normal to the beam 
axis. Invoking the Taylor hypothesis as in [1,2, 3, 8] we replace n(a,p') in (2.14) by 

n(a, p' - v(a)t) 


so that the response is now a function of time t as well. Thus (2.14) becomes: 
k 2 da 


[ x k* da f 


ex P nfg, p - v(a)t) \dp'\, (3.1) 


r) Jr 2 2 7 (a)(x - a) 

p € R 2 . 


Log Amplitude Scintillation 

The log-amplitude scintillation is given by: 

X (x,p,t) = log|£(x,p,t)| = log|E 0 (x)e^’ t) | 

= log|£o(x)| + R eip(x,p,t). (3.2) 

We are only interested in the part due to the wind: 

R eip(x, p , t) 

which we shall continue to denote 

P.0- 
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(3.3) 


Thus we have for the log amplitude scintillation due to wind: 

x(x, p, t) = f do [ h R (x - o,o,p ,p)n{o,p' - v{o)t) \dp'\ 
Jo Jr? 

i t 2 


where 


h r = Re 


2tt 7(a)(i - 


1 ik \p' - l(a)p\ 2 \ AS 

— eXP T 7 W(x-a)j' (3 ' 4) 


This is our basic “scintillation response” space-time field, which we shall use to deduce 
the spectrum and covariance functions. 


Scintillation Spectrum 

From (3.3), we see that for fixed x, and each p € R 2 , x(x,p,t) describes a tem- 
porally stationary Gaussian random process in — oo < t < oo. We proceed now to 
calculate the corresponding spectral density. For this purpose we first calculate the 
covariance function. 


R x (pup 2 ,t) = E[x(x,pi,t + s) x(x,p 2 ,s)] 


I J L da hR ^ X ~ a ' P> ' Pl ^ hR ^ X ~ a ' a ' p "' p2> > 

• R{p’ - p" - v{a)t) 


- L£ d ° (/*m*-w,*)^v) 

• ^h R ( I -c,a,p\p 2 )e- M ^'Up''y- M ^^‘Q(\)\dX\ (3.5) 


, h(x - a, op ', pi) + h(x - a, o, p', pQ 

o,o, p , pi) 2 

where 

, x k 2 1 ik \p' - x(o)pi,p' -l{o)pi\ 

Hx-a.a.p.n) = 


Now for Rep > 0: 

— / e^'^exp— | dp '| = exp(-27r 2 A 2 /x + 27ri[A 1 m]). (3.6) 

27T/X 2 p 
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f e 2ni[x,p " ] h R {x - a, a, p”, p 2 ) \dp " | 

JR* 

1 [ — 47t 2 A 2 17((t)(i — a ) , 7 r 

= - tfcexp exp-27rt[A, p 2 ]7( cr ) 

4 n 2 \ 2 i'y(o)(x - a) f , , , 

— ifcexp — exp — 27ri[A, p 2 ] 7 (a) . 

Zk 

Hence a little calculation leads to 

R x (pi,P2,t) = f X da f e- 2 Ki[x ' v{a]]t Q x {\,o,Pi,p2)QW |<*A| (3.7) 

Jo J pJ 
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where 


Q X (A, a, pi, pz) = j ^exp47r 2 A 2 7 / (a)^j^^exp27ri[A,p 1 - p 2 ]l R (°) 

■ [2 cosh 27r [A , pi + P2]l I i cr ) 

_ e ^ n2 ^ 2 ^f R ( a ) £ T Z e ^ 2 k[\ pi-p2]7 / ( t7 ) 

-47T 2 iA 2 7 (a) 2 tt[A, pi -P2]7 7 (&) 

— e « o i 

= exp ^47r 2 A 2 7 ; (a) T ^ — ^ exp(27ri[A, pi - P 2 ]l R (^)) 

2 cosh ^[A.px 4- P2]7/(^)) 

- 2 cosh ^ a i ~ - Shr[*. Pi - Pull/W) . 

0 < a < x. (3-8) 

Prom (3.7) we see that the scintillation intensity at p\ 

R x {p,p,0) = E[{x(x,p,t)) 2 } 
does not depend on the wind velocity. Defining 

Q*(A,a,p) = Q x (X,a,p,p) 

we have: 

Q x (X,a,p) = y^exp4?r 2 A 2 7 

2 cosh 47 t[A, p]7 / (ff) — 2cos47r 2 A 2 7 fl (a) — — — . (3.9) 

The cross-spectral density P{pu P2< f) is defined by the Fourier transform: 

P (pi i P2i / ) = / e~ 2mft R x (pi,p 2 ,t) dt. 

j—oo 

Thus the spectral density matrix of the 2 x 1 scintillation process at any two points 
Pit Pi'- 


Reip{x,pi,t) 
Reij)(x l p2, t) 


Pit P2 €■ R 
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for fixed x is given by 

P(pl,puf) P{p\,P2,f) 

P{p% Pl,f) P(p2,P2,f) 

We shall now proceed to calculate this matrix function, which is the main contri- 
bution of this paper. First we shall specialize to the plane wave case to compare with 
known results. 


Spectral Density: Plane Wave Case 

As we have indicated, the plane wave case is obtained by setting 


7 ,(<0 = 0; 7 *(*) = 7(<0 = !• 

Substituting into (3.8) we have: 


^.2 r 


(7, pi, P2) 2 

and in turn we have: 


1 — cos47r 2 A 2 


x — a 


exp27ri[A,pi - p 2 ] 


(3.10) 


jfc 2 


RAP uK.t) = I’ da J R e— -re 


27 ri[A f v(<7 )]t p 2m[\ } pi-p2] 

R 2 2 


[ 9 9 (x - a) 

• 1 — cos47t 2 A 2 - 


Q( A) \dx\ 


= Jo da f(-) (t) A>/ o(27rA|v(a)t + (p 2 - Pi)l) 


1 — cos47r 2 A 2 


(x-a) 


Q( A) dA. (3.11) 


This depends only on the difference in position (pi - p 2 ) and 

x(*>p. 0 


is thus both spatially and temporally stationary. This function is noticeably different 
from that for the case of the spherical wave derived by Lawrence et al. in [9]. 
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Specializing to p x — p% = p , we have: 


R x (p, P, t) = j (jj 1 - COs47T 2 A 2 X - - Q 

rX f k^\ f°° 

= y 0 da (2tt) y J Q \Jo{2n\\v(<j)t + Pi- P 2 \) 


Q( A) |dA| 


• 1 — COS 


47r 2 A 2 (x — a) 


Q( A) rfA. (3.12) 


Note that the covariance function does not depend on p; and depends only on the 
magnitude of the wind velocity | v(cr ) | . 

We shall denote corresponding spectral density P(f,p) by P(f): 


Then 


where for / > 0, 


R x (p,p,t) = r e M,, PU) if 

J — oc 

p (/) = f x rrv\ Pi ( a ' rAA) da ' 

Jo |w(a)| V K a )l / 


{aJ) = t T ~ cosi ' vL ir) Q{x) 


(3.13) 


(3.14) 


If the velocity is a constant: 

|v(a)| = v 

the spectral density formula (3.14) simplifies, since we can perform the integration 
with respect to a. In fact 

p(f) = - Pi (-) 


where 


pm - — H kQW (\ 

2 /f s/X‘ - p \ 


sin47r 2 A 2 

4,n= (i) 


(3.15) 


agreeing with the result obtained by Ishimaru [2], who also obtains therefrom the 
estimate: 

Pi(f) ~ / _8/3 for large/. 
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Cross-Spectral Density: Plane Wave 

Let us next consider the cross-spectral density. We have: 

P(f,Pi,p2) = [ P(<r,f,Pi,p2) 

Jo 

where 

P (<7, /, Pi . Pi) - / e~ 2wift R x {a, pi, p 2 , t ) 
is the “distributed” spectral density, where 

2x fc 2 

R x {o,Pi,p 2 ,t) = — — J Q XJ 0 {2n X\v(a)t + p\) 


dt 


^1 — cos4tt 2 A 2 — - — ^ Q(A) dX 


where 

P = Pi~ P2- 

By the Bessel function summation formula: 

J 0 (2nX\v(a)t + p\) = J 0 (27tA |v(a)| t) J 0 (27rA|p|) 

OO 

+ 2^2 J m (2n X |v(cr ) 1 1) J m (2n X\p\) cos(ro7r 4- <p(o)) t > 0 

l 

= Jq(2itX |v(<t)| t)Jo(2n X\p\) 

OO 

+ 2^>2 J m (2n X |u(a)|t)J m (27rA|p|)cosm<6(a), 

l 

where the angle is defined by 

M o),p] = |v(a)| |p|costf>(a). 

We see that the cross-spectral density depends now on the angle 4>{o), although we 
can write 

p (o,f,pi,p2 ) = p i | r ( a )|’ pl,p2 ) 

where the subscript “1” again stands for the case |v(ct)| = 1, and we note that 
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/ oo fOO 

J 0 ( 2ttA 1 t + p\)e~^ ift dt = / J 0 (2tt Xt)e~ 2mft dt J 0 (2* A|p|) 

-CX> J OO 

00 roo 

+ 4t j cos ft J-2m(2n\t) dt d2 m (27r A|p|) 

°° roo 

- 4 t ^ sin 27 r/t J 2 m-i( 27 rAt) dt d 2 m-i( 27 r A|p|). 

Now by [6, p. 405], for A, / > 0: 

f cos 27r/td2m(27r At) dt = 0, 


roo 

/ sin 27r/t J2m-i(27r At) dt = 0, 
do 


Hence 


+ 


= 0, 


A < / 

cos 

(2m arcsin (f)) 

A > / 


2l Ky/\*-f* ’ 

= 0 , 


A < / 

sin 1 

((2m — 1) arcsin (^^ 

A > /■ 



:J 0 ( 27 rA|p|)AQ 2 (A,a) dA 



1 * 00 cos (2m arcsin (^)) 


if vA* - /« 


■ d 2 m ( 27 rA|p|)AQ 2 (A, 0 ) dA 



.17) 


<* /oo sin ((2m — 1) arcsin ({ ) J 


^ 2 m-i( 2 jtA|p|)A< 3 2 (A,<T) d A 


(3.18) 


where 


Q 2 (A,a) = y 1 - cos ( 47r2 A 2 — £— )] Q(A). 

This result would appear to be new with this paper. We note that 

Pi(<r,f,Pi,Pi) = Pi(v, /.P 1 .P 2 ) • 


(3.19) 


19 



Spectral Density: General Beam Wave Case 

Let us now consider the general beam wave case. We shall calculate the spectral 
density of the scintillation at the point p. Let 

R(p,t) = R x (p , p,t). 


Then 

R(p,t) = f do R(o,p,t) do 
Jo 


where 

Since 

R(o,p,t) = / e- 2 ^ x ' v ^ t Q x (X,o,p)Q(X) \dX\. 

(3.20) 


Q x (X,o,p) > 0, 


we see that R(o,p t t) is a covariance function in t, — oo < t < oo. Also 

R(p, 0) = ( X do [ O x (X,o,p)Q(X)\dX\ (3.21) 

Jo Jr 2 

and as a function of p, depends only on |p|, and is a minimum at p — 0, since 

cosh 47r 2 [A, pYlj{o) > 1. 

As noted earlier (3.21) does not depend on the velocity v(-). Let <f>(a) denote the 
angle between the vector p and v(<r), so that: 

[v(a),p] = |v(ct)| \p\ cos 0(a). 


Then 

R(a,pj) = r° j*' e -2,ri l A ll 1, ( <T )l tcostf A<5 x(A, o ) 

• [cosh(47r 2 |A| |p| cos(0 — <t>(a))'y j(o)) 

— cos ^47r 2 A 2 7 R (a) — - — dX dO (3.22) 

where we have used the notation for short: 

Qi(A,a) = y Q(A)exp ^47r 2 A 2 7 / (a)^-^-^ . (3.23) 
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Now we recall the relation [6, p. 358] 

i r 27 T OO 

— jf e *«"'e bce *'-+) dO = J o (a)J o 0) + 2£ j k ( a ) J k (() cos kj. (3.24) 
Letting 

a = 27r | A| |u(<j)| t 
b = 4tt | A] |p| 7 / (cr) 


we have, using (3.24), that 


e -2ni\\\\ v ( a )\tc° a 9 cosh (47r 2 |A| |p| 7/ (a) cos(0 - tf(a))) dQ 

OO 

J 0 (a)I 0 (b) + 2^(-l) fc J2*(a)/2*(&)cos2fc0(<7) 


= 27T 


= 27rJo(a) 


for p — 0. 


(3.25) 

(3.26) 


Also 


Since 

we see that 


R(a, 0, t) = 47i - j Jo(27r|A| |n(a)|t)A<3 1 (A,a) 

• ^1 - cos ( 4 ^ 2 A 2 7 R (^)— j—)) d<J - 

Q x {X,a,p) - Q x (X,a,0) > 0, 


(3.27) 


R(a,p,t)-R(a,0,t) = / (Q*(A, a, p) - Q x (X, a, 0)) |dA| 

J H '? 

is a covariance function in f, — oo < t < oo; and further, using (3.25), (3.26), 


= 2t r jT AQ^A.a) 


d X. (3.28) 


Jo(o)(fo(&) — 1) + 2 YX-l) k J*(a)I*{b) cos2k<t> 

l 

Here we may decompose R(cr,p,t) as: 

R(a,p,t) = R{cr,0,t) + R(a,p,t) 

where R(o, p , t) is given by (3.28). Correspondingly we decompose the spectral density 

P{f,p)=f P(a,f,p)da, -oo < / < oo, 

Jo 
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setting 


P{a,f,p) = P(<r,f,0) + P(a,f,p) 
P(cr,f, 0) = 2 jf R(a,0,t) cos2n ft di 

ro o _ 

P{o,f,p) = 2/ R{a,p,t) cos 2nft dl 
do 


Now 


P(<7,/ ' 0) Pl ("’ |®(ff)| 

where the subscript “ 1 ” corresponds to the special case 

|v(a)| = 1. 


■°) 


Using (3.16) we have for 0 < /: 


Afa/.O) = 


r° AQi(A.a) 

J f 


COS 


(4tt 2 A 2 7 r (ct) : 


Similarly we have 


P(<r,f,p) 



where 


/•OO 

P\{o,f,p) = 

J J 


AQi(A, a) 


where 


(/o(ft)-l) + 2^(-l) A: cos 
1 


^2fc arcsin 



b = 47r 2 A|p|7 / (g). 


Hence finally we have: 


P(v,f,p) 



(3.29) 



hk(b) cos 2k<p(a)\ d A 
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and for / > 0: 


P\(oJ,p) 


For p = 0: 


f°° \Qi(X, a) 

Jf 


1 0 (i>) + 2y^(-l) fc cos ^2fcarcsin T&W cos 2k<t>(a) 


~ cos (47r 2 A 2 7 *(a) 


d X. 


(3.30) 


r['- 


XQ(X) 


oo 

f 2 ' T^T 2 


COS 


' 47r 2 A 2 7 B (a)(a: - a) 

k 


1 


exp 


/4tt 2 A 2 


V 


7 i(a){x ~ <0^ 


da (3.31) 


If we specialize to the plane wave case, 

Pi(<b/,p) = Fi(a, /, 0) 

and in particular (3.31) simplifies of course to (3.15). Put another way, (3.31) is the 
generalization of (3.15) to the beam wave case. 


Spectral Moments/Path Weighting Function 

As we have noted, the primary impetus in studying scintillation was in “remote 
sensing,” or estimating, the wind velocity. Thus Lee and Harp [8] noted that the time 
derivative at time zero of the space time cross-covariance function could be expressed 
as a linear function of the wind velocity. They were working with the spherical wave 
model but we can see this as well from (3.11) derived for the plane wave. Thus we 
can calculate that 

X *xte.P2,l)l„o = [\v(o)\W(o,p)io (3.11a) 

at JO 

where 

W{a, p) = 47T 2 A 2 Ji(27tA|p|) ^1 - cos4tt 2 A 2 — - — ^ Q(A) dX 
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and note that it does not depend on the wind velocity. Thus it serves as a “path 
weighting function.” 

One disadvantage of this technique involves taking the slope of the cross correlation 
function which is an off-line operation, and further the path-weighting function is 
equal to zero when p = 0. 

We shall show that it is possible instead to work with the spectral moments. 

Given any spectral density P(-) (of any real-valued process), we may define the 
“spectral moment of order n,” by 


75 = jo rrpiM. 
1 Jo°°P(f)df 


(3.32) 


yielding a measure of the spectral spread. We can express even-order spectral mo- 
ments in terms of derivatives of R{-), the corresponding covariance function. Thus 
for n = 2, we have 

_ r?"fm 

(3.33) 


**7* = ~ R " i0) 


fl(0) 

On the other hand we should note that all odd derivatives of R(-) vanish at the origin. 
Here we shall calculate all the spectral moments directly from (3.30). Thus we 

Jo M^)P a n(<bP) da 


f n (p) = 


Jo fl o(<b P ) d<J 


(3.34) 


where 


roc 

a„ = / rPi(<J,f,p)df 

Jo 

= ^T /2 (sin 9) n do'j 


■ J o A n+1 Q!(A,a) (/ 0 (6) - cos ^47r 2 A 2 7 ft (a)^-^^ dX 


r ir/2 


) 


+ 2^(-l) fc / (sin0) n cos2 k9 dO 


roo 

/ A n+1 Qi(A, a) (hk{b) cos2k<f>(a)) d A. 
Jo 


(3.35) 
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In particular: 

a 0 (o,p) = l£°\QiM(lo(b) - cos(47r 2 A 2 7 R (a)^)) dA 

ai(<r,p) = J q \ 2 Qi(\,<t) (l 0 (b) - cos ^47 t 2 A 2 7 r ((t) 

a 2 (a,p) = jJ q A 3 Q!(A,a) ^/ 0 (6) + / 2 (*>) cos2<£(a) 

- cos ^47r 2 A 2 7 H (g) 1 ^ dX. 

We may express (3.32) as 

/ n (p) = / |t;(a)| n W n (g,p)d(7 

do 

where W n (<7, p) is the “path weighting function,” and is given by 

a n (a,p) 


W n (a,p) = 


J^a 0 (a,p) da 


(3.36) 


(3.37) 


We note that for the general beam wave case the path weighting functions de- 
pend on the angle <p(a) for nonzero p. However this dependence disappears when we 
specialize to the case of the plane wave: 

a n (a,p) = y (j ' (sinfl)" d0^ jf A n+1 Q(A) ^1 - cos4tt 2 A 2 —J— ) <*A (3.38) 

and in particular: 


f a 0 (a, p) 
Jo 


da = 


7r fc 2 r°° 


J°° AQ(A) (/ (l — cos47r 2 A 2 |) da^ dA 
/•°° / sin47r 2 A 2 f\ 

L xm [ l - -wfj " 


4 

7rA: 2 x P 00 


Thus the path weighting functions do not depend on p, for the plane wave case. 


Zero Crossing Rate 

The spectral moment of order two is of particular importance because of Rice’s 
theorem [18] relating it to the zero-crossing rate of the process: 


Average zero crossing rate = 2 . 


(3.39) 
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Thus yields, in our notation: 


No{p ) 2 = 4 f |v(cr)| 2 VT 2 (a,p)da (3.40) 

Jo 

where No(p) is the average zero-crossing rate at the position p. The zero-crossing 
rate can be readily measured on-line by analog instrumentation. Thus (3.40) offers a 
postential alternate technique for measuring wind velocity. 


4. Random Wind Model 

In this section we model the wind velocity vector field v(a), 0 < a < i, as random 
— and more specifically, as Gaussian distributed with mean m(a) and nonsingular 
covariance matrix R(a). We seek to calculate the average spectral density: 

P(f,p) = E[P(f,p) ]■ 


Rather than use (3.30) we calculate 


E[R(p,t)\ 


where by (3.20) 


R(p,t) = [da [e- 2 ^ x ' v ^ t Q x (X,a,p)Q(X) |dA|. 


Hence 


E[R(p,t)} = [da [E[e- 2 ^ x ' v ^ t }Q x (X,a 1 p)Q(X) |dA| 

and we see that only the first-order distributions of v(-) are involved. Now 
£[ e -27 r i[A,t>(<7)|tj _ exp(-27r 2 [R(a)A, A]t 2 + 27ri[A,m(cr)]0 

and hence 

r e ^nE{e~^ i[XMa)]it ] dt 

7 — 00 


v^ry[fl(»)A,Ai exp 2 


-1 (/ - [A, m(a)])‘ 


(4.1) 

(4.2) 


(4.3) 


(4.4) 
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Hence 


P(f,p) = f X da [ —4 = Q x (X,a,p)Q(X) 

Jo J* ^y/[R{a)X, A] 


-1 (/ - [A, m(a)]) 2 IJA( 

• eXP T [R(a)X— m 


(4.5) 


Diagonal Covariance 

While the mean m(a) can be arbitrary, the covariance matrix model may be 
reasonably simplified to be diagonal: 

R(a) = d 2 (cr)/ 2x2 . ( 4 - 6 ) 


In this case we can derive a useful alternate form for (4.5). Thus 
E[R(p,t)\ = J*da Q X (A, a, p)Q(X) 

• exp[27ri[A,m(a)]t - 27r 2 A 2 d 2 (a)t 2 ] |dA| 

= 27 r J da j AQi(A, a) (exp — 27r 2 A 2 d 2 (a)t 2 ) 

[j 0 (27r A|m(a')| |t|)/o(4TrA|p|7 7 (a)) 

OO 

+ 2^(-l) fc 7 2 jk(27rA|m(a)| \t\) 

1 

• / 2 /k(47r A|p| 7 7 (a)) cos 2k<p(a ) 


Q O /v 

— cos47r A 1 r {o) 


x — a 


d A, 


where 4>{a) is the angle between p and m(a). Hence, for |m(a)| ^ 0, 
J e _27r * /f (exp-27r 2 d 2 (a)A 2 l 2 ) J 2 k{ 27rA|m(a)| |«|) dt 


1 1 rA|m(<r)| 1 


= ~-f 

2n X Jo 


V^d(a) 


( -1 U-uf , -1 (j±»l 

K eXP 2 A 2 d 2 (a) eXP 2 A 2 d 2 (a) 
cos (2k arcsin fa)) 


\j A 2 |m(cr)| 2 - J/ 2 


du. 


(4.7) 
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Hence 


P(f,p) = f da J d\ ■ XQx(X,a) 

■ [/ 0 (47rA|p|7 / (a)) 4- 2 hk{^ Mp\l ^)) 


• cos ( 2k arcsin 


A|m(a)|)) 


cos 


2 k<t>{a) 


- cos47r 2 A^7 R {a) 

•|m(<r)|A \ 1 


X — a 


\Phi A d{a) 


-i (/ - v ) 2 

exp -^v 7 + exp 

i 


-i (/ + 1^) 2 \ 

2 A 2 d 2 {a) ) 


A 2 |m(a)| 2 - i/ 5 

We note that as d{a) -* 0, the last factor in parentheses: 


dv . 


(4.8) 


— + — ==:==^= for / < A|m(a)| 

\J A 2 |m(cr)| 2 - / 2 

= 0 for / > A|m(cr)| 

and thus (4.8) agrees with (3.30) for d(a) = 0. We may interpret the deterministic 
case as mean wind of variance zero. 


Case of Zero Mean Wind 

On the other hand, one advantage of the random wind model is that we can 
consider the case where the mean wind is zero, unlike the deterministic case. Thus 
let 

m{a) = 0, 0 < a < x 

and further consider a coaxial detector: 


P = o. 
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Then (4.5) yields: 


P(/,°) = l* da 


d(a) 


■ exp 



d X. 


(4.9) 


This has no analog in the deterministic case. It is interesting to note that the spectral 
density is now a monotone decreasing function of the frequency. We can simplify 
further if we take wind-covariance to be constant: 


d(a) = d, 0 < a < x. 

Specializing to the plane wave case, (4.9) simplifies to: 


— Jfc 2 X f°° n/27T /- 1 / 2 > 

p(/ ' 0) = — Jo — exp {TM, 


l - 


in 2 A 2 (f) 

4„n= (f) 


) 


Q( A) dX. (4.10) 


If further as in [9] we invoke the “inertial subrange approximation” for the Kolmogorov 
spectral density: 

Q(A) ’ 1 


( J , + 4,^) 11/6 ( 2 ^) U/3 

we can simplify (4.10) further to: 

_ k 2 x f°° 1 /— / 2 \ 

P(/ ' 0) K — Jo — (2^73 exp (^j 

and evaluating the integral on the right yields: 


d X 


j ) ' <4 - n) 

Note the similarity to the estimate in the deterministic case (cf. (3.15)). 


Spectral Moments/Path Weighting Functions 

We shall consider only the most useful case: the spectral moment for n = 2 and 
the corresponding weighting function. Thus using (4.2) we have: 
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2 iT /2jP(/,p)d/ 

= f da [ E([X,v {a)} 2 )Q x (X, a, p)Q(X) \dX\ 

Jo Jr 2 

= J* da J Ri (\[X,m(a )}\ 2 + [R(a)X,X]}Q x (K(r,p)Q(X)\dX\. 


For the term containing m(a) we have only to replace v(a) by m(a) in (3.33) for 
n = 2. For the second term, using (3.24), we have that 

[ jR(a)X,X]Q x (X,a } p)Q(X) |dA| 

J R? 

= 2nd 2 (a) I A 3 Qi(A, a) J/ 0 (47rA|p|7 / (a)) - cos ^47r 2 A 2 7 H (a)^-^-^ . 

roc 

2/ P(f,p)df = E[R(p, 0)] 

Jo 

and by (4.7) yields: 

rX roo 

E[R{p, 0)] = 2tt J da J XQi(X, a) 


Now 


I 0 (b) - cos ^4n 2 X 2 y R (a)^Y~^ dX - 


Hence 


757-7 f£°f 2 P(f.p)Jf 

/W JS°P(f,P)df 


can be expressed as 
where 


= flip) + !S(p) 

TU7) = f \m(o)\ 2 W 2 {o,p) da 

Jo 

where W 2 (a, p) is given by (3.36). And 

fo(p) = j (d 2 (o)) da W 2 {a,p) da. 
In other words the spectral moment of order two: 

IHp) 
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for the random wind is the sum of two terms, the first corresponding to a deterministic 
wind velocity equal to the mean wind velocity and the second corresponding to a 
determinstic wind velocity equal to the “sigma 11 of the random wind: 

d(a) = y/d 2 (a) . 

The path weighting function of course does not depend on the statistics of the wind 
velocity. Note that the variance term represents an additive error term. This “addi- 
tive” property does not extend to higher order spectral moments. 
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